How do I solve a differential equation in Mathematica?

Note: This is not a primer on differential calculus, differential equation solving methods or numerical accuracy. If you are interested in any of these topics, please consult the appropriate literature. Some good books on these subjects can be found in the Mathematica bookstore.

One important advantage of Mathematica is the ease with which users can solve many families of differential equations either numerically or symbolically.  for example.

In many cases, it is a good idea to try the symbolic solver DSolve first. With the constant improvements made to Mathematica functions, DSolve can often find a complete analytic solution, even for higher-order or non-linear differential equations.

Symbolic solutions

What kinds of problems can DSolve solve?

DSolve can solve linear ordinary differential equations of any order with constant coefficients. It can solve also many linear equations up to second order with non-constant coefficients. DSolve includes general procedures for many of the nonlinear ordinary differential equations whose solutions are given in standard reference books such as Kamke. It can also can find general solutions for linear and weakly nonlinear partial differential equations. Truly nonlinear partial differential equations usually admit no general solutions.

Ordinary differential equations

All most users will have to do is just type in the differential equation or the system of differential equations. As opposed to many numerical differential equation solvers, Mathematica does not require users to transform higher order equations to systems of first-order equations by hand. Below are some small examples. For more information on DSolve, see the documentation, especially the Further Examples.

Ordinary differential equations without initial conditions

In this case, you can usually just type in the equation and evaluate the DSolve statement.

DSolve[y'[x]-x y[x]==0,y[ x],x]

DSolve also works for higher order equations,

DSolve[y''[x]-x y[x]==0,y[x],x]

DSolve[Derivative[3][y][x]-x y[x]==0,y[x],x]

DSolve[Derivative[4][y][x]-x y[x]==0,y[x],x]

DSolve[Derivative[4][y][x]-x (y'[x])==0,y[x],x]

as well as for first order equations with non-linear coefficients.

DSolve[(y^( (1)))[x]-(x^2) ((( y[x]))^2)==3 (x^3),y[x],x]

And just for the heck of it, let's give Mathematica this differential equation:

DSolve[(y^( (41)))[x]-a x (y^( ′))[x]==0,y[x],x]

Ordinary differential equations with initial conditions

These equations are handled pretty much the same as ordinary differential equations without initial conditions. Just enter the initial conditions as part of the list of equations.

DSolve[{Derivative[1][y][x] - x*y[x] == 0, y[0] == 0.1}, 
  y[x], x]

Systems of ordinary differential equations

Solving systems of ordinary differential equations is the same as solving a single ode. Simply enter the list of equations and run DSolve on it.

DSolve[{y[x]==-z'[x],z[x]==- a y'[x],y[0] == 1,z[0] == 0},{y[x], z[x]},x]

Partial differential equations

The theory of partial differential equations is less fully developed than the one of ordinary differential equations: specifically, there is no equivalent of the Cauchy-Kovalevskaya theorem. This means that Mathematica can only solve a subset of differential equations symbolically. Specifically, DSolve can find general solutions for linear and weakly nonlinear partial differential equations. Most nonlinear partial differential equations do not allow for general solutions. In these cases it is advisable to use the tools in the Calculus`DSolveIntegrals` standard package to attempt to find complete integrals or differential invariants.

DSolve can find general solutions to partial differential equations like the wave equation in the following example. Once again,  just enter the equation and call DSolve.

DSolve[-Derivative[0, 2][y][x, t]+(c^2) Derivative[2, 0][y][x, t]==0,y[x, t],{x, t}]

Numerical solutions

In Mathematica, the NDSolve command works consistently for first order ordinary differential equations, partial differential equations, systems of differential equations and higher order differential equations. Additionally, it's default method, Automatic, automatically switches between Gear and Adams, depending on stiffness. Another possibility for equations which are not stiff is RungeKutta. For some problems, the Runge-Kutta method can find the solution faster than the other methods.

Higher order differential equations

NDSolve requires the differential equation(s), a sufficient set of initial conditions, and the range to solve for as inputs. The /. shortcut in the Plot command inserts the solution for the function y(x).

solution=NDSolve[{Derivative[3][y][x]+y''[x]+y'[x]==-y[x]^3,y[0] == 1,y'[0]==y''[0]==0},y,{x, 0, 20}];Plot[y[x] /. solution, {x, 0, 20}];

Boundary value problems

NDSolve can solve a significant number of boundary value problems. Just add the boundary conditions as additional equations.

L = 10;α = 1.1;solution=NDSolve[{(u')'[x]+(α^2) u[x]==0,u'[0]+I u[0]==0,u'[L]-I u[L]==-2 I (E^(-I L))},u,{x, 0, L}];Plot[Evaluate[Abs[u[x] /. solution]], {x, 0, L}];

Systems of differential equations

In this example, Mathematica solves the equations of a Lorenz attractorin one step. In other packages, the user would have to do significant setup work.

solution = 
        NDSolve[{Derivative[1][x][
                t] == -3*(x[t] - y[t]), Derivative[1][
                  y][t] == 
              26.5*x[t] - y[t] - x[t]*z[t], Derivative[1][z][
                t] == x[t]*y[t] - z[t], x[0] == z[0] == 0, 
            y[0] == 2}, {x, y, z}, {t, 0, 17}];

Here is a graph of the solution.

ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. solution], {t, 0, 17}, 
      PlotPoints -> 1000]

Partial differential equations

In this example, NDSolve solves the heat equation with the left end (x = 0) held at fixed temperature 0, the right end (x = 1) held at fixed temperature 1, and an initial heat profile given by a quadratic in x.

solution=NDSolve[{(∂_t) u[x, t] == (∂_(x, x))u[x, t], 
            u[x, 0] == 6 x - 6 (x^4) + (x^2), u[0, t] == 0, u[1, t] == 1}, 
          u, {x, 0, 1}, {t, 0,  .5}];
Plot3D[
        Evaluate[u[x, t] /. First[solution]], {x, 0, 1}, {t, 0, 
          0.5}];

A more interesting example is when time-dependent boundary conditions are added. For example, entering the following produces a plot of a solution with the temperature at the left edge varying sinusoidally.

solution=u/.First[NDSolve[{\(∂\_t u[x, 
                          t] == \(2\ ∂\_{x, 2}u[x, t]\)\/\(9\ π\^2\)\),u[x, 0] == 0,u[0, t] == Sin[t],Derivative[1, 0][u][1, t]==0},u,{x, 0, 1},{t, 0, 6 π}]];Plot3D[solution[x, t], {x, 0, 1}, {t, 0, 6 π}, 
        PlotRange -> All, PlotPoints -> 25];

Systems of partial differential equations

NDSolve can also handle systems of partial differential equations. The notation {u, v}/. just means "replace the functions u and v with the results of NDSolve".

solution={u, v}/.First[NDSolve[{(∂_t) u[t, 
                          x] == (u[t, x] - 
                              1) (-(16 (v[t, x] - 1)) - 
                              2 t + 
                              16 t x) + (∂_x)((v[t, x] - 
                                  1) (∂_x) u[t, x]) + 
                        10 x (E^(-4 x)),\(∂\_t v[t, x] == 
                      x\^2 - 2\ t + ∂\_x u[t, 
                            x] + ∂\_{x, 2}v[t, x] - 
                        10\ t\ E\^\(\(-4\)\ x\) + 4\ u[t, x] - 
                        4\),u[0, x] == 1,v[0, x] == 1,u[t, 0] == 1,v[t, 0] == 1,3 u[t, 1]+Derivative[0, 1][u][t, 1]==3,5 Derivative[0, 1][v][t, 1]==(E^4) (u[t, 1] - 1)},{u, v},{x, 0, 1},{t, 0, 2}]];

plot1  = 
        ContourPlot[solution[[1]][x, 
            t], {x, 0, 1}, {t, 0, 2}, ContourLines -> False, Contours -> 150, 
          ColorFunction -> Hue,  
          PlotLabel -> 
            StyleForm["u(t,x)", FontFamily -> "Times", 
              FontSize -> "18", FontSlant -> "Italic"], 
          DisplayFunction -> Identity];

plot2  = 
        ContourPlot[solution[[2]][x, 
            t], {x, 0, 1}, {t, 0, 2}, ContourLines -> False, Contours -> 150, 
          ColorFunction -> Hue, 
          PlotLabel -> 
            StyleForm["v(t,x)", FontFamily -> "Times", 
              FontSize -> "18", FontSlant -> "Italic"], 
          DisplayFunction -> Identity];

Show[GraphicsArray[{plot1, plot2}], 
      PlotLabel -> 
        StyleForm["Solutions", FontFamily -> "Times", 
          FontSize -> "24", FontSlant -> "Italic"], 
      DisplayFunction -> $DisplayFunction]

Checking results

The easiest way to check results in Mathematica is to substitute the solution into the original equation. This technique works for both symbolic and numerical solutions.

symsol=DSolve[{y'[x] ==  (1/(2 x y[x])), y[ .01] ==  .1}, y, x]

/.symsol replaces y[x] with the solution found by DSolve.

FullSimplify[(y'[x] ==  (1/(2  x y[x]))) /.  
        symsol]

The solution satisfies the differential equation.

numsol = NDSolve[{y'[x] ==  (1/(2  y[x])), y[ .01] ==  .1}, 
          y, {x,  .01, 1}];

solution1  = 
        Plot[Evaluate[(y'[x] -  (1/(2  y[x])) ) /.  
              numsol], {x,  .01, 
            1}, 
PlotRange -> {-(10^(-3)), 
              (10^(-3))}];

The solution found by NDSolve satisfies the equation in very good approximation. The rather high difference at the beginning is an artifact of the numerical differentiation. See the next section on how to even more accurate solutions.

Higher precision settings

One way to get a very precise solution of an ODE is to give a sufficiently high value for the WorkingPrecision option. Note that AccuracyGoal and PrecisionGoal default to 10 less than the value of WorkingPrecision when it is greater than $MachinePrecision.

Let's compare the result of an earlier calculation to the same one running with higher precision.

lowsol = NDSolve[{y'[x] ==  (1/(2  y[x])), y[ .01] ==  .1}, 
          y, {x,  .01, 1}];

highsol = NDSolve[{y'[x] ==  (1/(2  y[x])), y[ .01] ==  .1}, 
          y, {x,  .01, 1}, WorkingPrecision -> 22];

The following plot shows the logarithm of the difference between the numerical solution and the differential equation with the default setting of WorkingPrecision (black) and with WorkingPrecision set to 22 (red).

Plot[
        Evaluate[-Log[
              Abs[{(y'[x] -  (1/(2  y[x])) ) /.  
                    lowsol, (y'[x] -  (1/(2  y[x])) ) /.  
                    highsol}], 10]], {x,  .01, 1}, 
        PlotStyle -> {GrayLevel[0], RGBColor[1, 0, 0]}, 
        AxesOrigin -> {0, 0}];

As you can see, the numerical errors are significantly smaller when using the higher working precision, the rather high first value is an artifact.

Other places to look

Calculus`DSolveIntegrals`

Most nonlinear partial differential equations do not allow for general solutions. In these cases it is advisable to use the tools in the Calculus`DSolveIntegrals` standard add-on package to attempt to find  complete integrals or differential invariants.

Integral Transforms

Mathematica includes many functions for calculating integral transforms. Especially Fourier and Laplace transforms are powerful tools for finding analytical solutions to differential equations. See the sections about Fourier and Laplace transforms in The Mathematica Book for more information.

Lie Symmetries

Gerd Baumann of the University of Ulm, Germany, has submitted a number of packages that deal with Lie and Lie-Baecklund symmetries. Lie symmetry methods can be used to construct solutions for linear and especially for nonlinear systems of differential equations.